% gets the result allowing precautionary savings
    
function [C_PS,A_PS,CE_PS,C_PS_PHI,A_PS_PHI,CE_PS_PHI] = analysis_PS(Wi,CAPH,MU,Ppaid_N,Pspot_N,iX,Nstates,rho,r,gamma,T,amax_phi,amax_st)    

    % Precautionary Savings with ST contracts
    [C_PS, A_PS] = C_Savings(Wi(:,1),MU,CAPH,gamma,iX,rho,r,amax_st);
    % here check that the max of A is lower than the amax_phi, but close
    max(max(A_PS))
    if max(max(A_PS))==amax_st*1000
        fprintf('increase max assets for ST');
    elseif max(max(A_PS))<0.8*amax_st*1000
        fprintf('decrease max assets for ST');
    end
    

    u = -1/gamma*exp(-gamma*C_PS);
    Num = mean(sum(u.*(iX~=(Nstates+1)).*power(rho,(1:T)-1)',1));
    Den = mean(sum((iX~=(Nstates+1)).*power(rho,(1:T)-1)',1));
    CE_PS= -log(-gamma*Num/Den)/gamma;

    %Precautionary Savings with PHI contracts
    [C_PS_PHI,A_PS_PHI] = C_Savings_PHI(Wi(:,1),Nstates,T,CAPH,Pspot_N,Ppaid_N,gamma,rho,r,amax_phi);
    % here check that the max of A is lower than the amax_phi, but close
    max(max(A_PS_PHI))
    if max(max(A_PS_PHI))==amax_phi*1000
        fprintf('increase max assets for GLTHI');
    elseif max(max(A_PS_PHI))<0.8*amax_phi*1000
        fprintf('decrease max assets for GLTHI');
    end    
    u = -1/gamma*exp(-gamma*C_PS_PHI);
    Num = mean(sum(u.*(iX~=(Nstates+1)).*power(rho,(1:T)-1)',1));
    Den = mean(sum((iX~=(Nstates+1)).*power(rho,(1:T)-1)',1));
    CE_PS_PHI= -log(-gamma*Num/Den)/gamma;
   
    
    end
    